(* Calculate and plot the Cdk1 nullcline *)

Clear[bkgd1,bkgd2,kcdc25,kwee1,ncdc25,nwee1,bigk1,bigk2,a,b,ec50,n,cdc25,wee1,ksynth,kdeg,cdk1star]

Solve[0==kcdc25*cdc25*(cyc-cdk1star)-kwee1*wee1*cdk1star+ksynth-kdeg*cdk1star,cyc]

cdc25=bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))

wee1=bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+(cdk1star^nwee1)))

kdeg=a+(((cdk1star^n)b)/((cdk1star^n)+(ec50^n)))

cyc/.%3

posfeedbackspeed=0.4;
kwee1=1*posfeedbackspeed;
bkgd1=0.20;
bkgd2=0.20;
kcdc25=2*posfeedbackspeed;
nwee1=3.5;
ncdc25=11;
bigk1=35;
bigk2=30;
n=17;
ec50=32;
synthdegspeed=.4;
ksynth=2.5*synthdegspeed;
a=0.025*synthdegspeed;
b=0.1*synthdegspeed;

ParametricPlot[{(cdk1star (a+((b (cdk1star^n))/((cdk1star^n)+(ec50^n))))+cdk1star (bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25-ksynth+cdk1star (bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+(cdk1star^nwee1)))) kwee1)/((bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25),cdk1star},{cdk1star,0,200},Frame->True,PlotStyle->{{Darker[Blue],Thick}},LabelStyle->{FontFamily->"Arial",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin, nM"],Style["Active CDK1, nM"]},PlotRange->{{0,102},{0,102}}]

(* Calculate and plot the cyc nullcline *)

ParametricPlot[{(ksynth/( a+(((cdk1star^n)b)/((cdk1star^n)+(ec50^n))))),cdk1star},{cdk1star,0,100},Frame->True,PlotStyle->{{Red,Thick}},LabelStyle->{FontFamily->"Helvetica",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin (nM)"],Style["Active cyclin-Cdk1 (nM)"]},PlotRange->{{0,102},{0,102}}]

Show[%23,%25]

FindRoot[(ksynth/(a+(((cdk1star^n) b)/((cdk1star^n)+(ec50^n)))))-(cdk1star (a+((b (cdk1star^n))/((cdk1star^n)+(ec50^n))))+cdk1star (bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25-ksynth+cdk1star (bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+(cdk1star^nwee1)))) kwee1)/((bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25),{cdk1star,30}]

y=cdk1star/.%27

x=(ksynth/(a+(((cdk1star^n) b)/((cdk1star^n)+(ec50^n)))))/.%27

ListPlot[{{x,y}},Frame->True,LabelStyle->{FontFamily->"Arial",FontSize->18},PlotStyle->{AbsolutePointSize[10],Black},PlotRangeClipping->False,PlotRange->{{0,100},{0,100}},AspectRatio->1,FrameLabel->{Style["Total cyclin, nM"],Style["Active cyclin-Cdk1, nM"]}]

ListPlot[{{x,y}},Frame->True,LabelStyle->{FontFamily->"Arial",FontSize->18},PlotStyle->{AbsolutePointSize[7.6],White},PlotRangeClipping->False,PlotRange->{{0,100},{0,100}},AspectRatio->1,FrameLabel->{Style["Total cyclin, nM"],Style["Active cyclin-Cdk1, nM"]}]

(* Here are the nullclines in the phase plane with the unstable fixed point plotted as well *)

Show[%26,%30,%31]

(* Numerically solve the two ODEs *)

tmax=500;

sol=NDSolve[{
cyc'[t]== ksynth-(a+b*cdk1[t]^n/(ec50^n+cdk1[t]^n))*cyc[t],
cdk1'[t]== (bkgd1+((((cdk1[t])^ncdc25) )/((bigk1^ncdc25)+((cdk1[t])^ncdc25))))kcdc25*(cyc[t]-cdk1[t])-(bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+((cdk1[t])^nwee1)))) kwee1*cdk1[t]+ksynth-(a+b*cdk1[t]^n/(ec50^n+cdk1[t]^n))*cdk1[t],
cyc[0]==0,
cdk1[0]==0},
{cyc,cdk1},
{t,0,tmax}]

(* Plot the solution as a time course *)

Plot[{cdk1[t]/.sol,cyc[t]/.sol},{t,0,500},Frame->True,PlotStyle->{Darker[Blue],Red},LabelStyle->{FontFamily->"Arial",FontSize->18},FrameLabel->{Style["Time"],Style["Active Cdk1 (nM)"]},PlotRange->{0,102},AspectRatio->1]

(* Plot the solution in the phase plane *)

ParametricPlot[{cyc[t],cdk1[t]}/.sol,{t,200,288},Frame->True,PlotStyle->{Black},LabelStyle->{FontFamily->"Helvetica",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin (nM)"],Style["Active cyclin-Cdk1 (nM)"]},PlotRange->{{0,102},{0,102}}]

Show[%32,%39]

(* And now the plot the trajectories based on DDEs *)

tau=15;

solDDE=NDSolve[{
cyc'[t]== ksynth-(a+b*cdk1[t-tau]^n/(ec50^n+cdk1[t-tau]^n))*cyc[t],
cdk1'[t]== (bkgd1+((((cdk1[t])^ncdc25) )/((bigk1^ncdc25)+((cdk1[t])^ncdc25))))kcdc25*(cyc[t]-cdk1[t])-(bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+((cdk1[t])^nwee1)))) kwee1*cdk1[t]+ksynth-(a+b*cdk1[t-tau]^n/(ec50^n+cdk1[t-tau]^n))*cdk1[t],
cyc[t/;t≤0]==0,
cdk1[t/;t≤0]==0},
{cyc,cdk1},
{t,0,tmax}]

Plot[{cdk1[t]/.solDDE,cyc[t]/.solDDE},{t,0,500},Frame->True,PlotStyle->{Darker[Blue],Red},LabelStyle->{FontFamily->"Arial",FontSize->18},FrameLabel->{Style["Time"],Style["Active Cdk1, nM"]},PlotRange->{0,102},AspectRatio->1]

ParametricPlot[{cyc[t],cdk1[t]}/.solDDE,{t,200,339},Frame->True,PlotStyle->{Blue},LabelStyle->{FontFamily->"Helvetica",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin (nM)"],Style["Active cyclin-Cdk1 (nM)"]},PlotRange->{{0,102},{0,102}}]

Show[%40,%45]

(* Now eliminate the ultrasensitive APC response *)

posfeedbackspeed=0.4;
kwee1=1*posfeedbackspeed;
bkgd1=0.20;
bkgd2=0.20;
kcdc25=2*posfeedbackspeed;
nwee1=3.5;
ncdc25=11;
bigk1=35;
bigk2=30;
n=1;
ec50=120;
synthdegspeed=.4;
ksynth=2.5*synthdegspeed;
a=0.025*synthdegspeed;
b=0.1*synthdegspeed;

ParametricPlot[{(cdk1star (a+((b (cdk1star^n))/((cdk1star^n)+(ec50^n))))+cdk1star (bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25-ksynth+cdk1star (bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+(cdk1star^nwee1)))) kwee1)/((bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25),cdk1star},{cdk1star,0,200},Frame->True,PlotStyle->{{Darker[Blue],Thick}},LabelStyle->{FontFamily->"Arial",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin, nM"],Style["Active CDK1, nM"]},PlotRange->{{0,102},{0,102}}]

ParametricPlot[{(ksynth/( a+(((cdk1star^n)b)/((cdk1star^n)+(ec50^n))))),cdk1star},{cdk1star,0,100},Frame->True,PlotStyle->{{Red,Thick}},LabelStyle->{FontFamily->"Helvetica",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin (nM)"],Style["Active cyclin-Cdk1 (nM)"]},PlotRange->{{0,102},{0,102}}]

Show[%63,%64]

FindRoot[(ksynth/(a+(((cdk1star^n) b)/((cdk1star^n)+(ec50^n)))))-(cdk1star (a+((b (cdk1star^n))/((cdk1star^n)+(ec50^n))))+cdk1star (bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25-ksynth+cdk1star (bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+(cdk1star^nwee1)))) kwee1)/((bkgd1+((cdk1star^ncdc25)/((bigk1^ncdc25)+(cdk1star^ncdc25)))) kcdc25),{cdk1star,30}]

y=cdk1star/.%

x=(ksynth/(a+(((cdk1star^n) b)/((cdk1star^n)+(ec50^n)))))/.%%

ListPlot[{{x,y}},Frame->True,LabelStyle->{FontFamily->"Arial",FontSize->18},PlotStyle->{AbsolutePointSize[10],Black},PlotRangeClipping->False,PlotRange->{{0,100},{0,100}},AspectRatio->1,FrameLabel->{Style["Total cyclin, nM"],Style["Active cyclin-Cdk1, nM"]}]

ListPlot[{{x,y}},Frame->True,LabelStyle->{FontFamily->"Arial",FontSize->18},PlotStyle->{AbsolutePointSize[7.6],White},PlotRangeClipping->False,PlotRange->{{0,100},{0,100}},AspectRatio->1,FrameLabel->{Style["Total cyclin, nM"],Style["Active cyclin-Cdk1, nM"]}]

Show[%69,%70]

(* We've done the same for the two stable fixed points and plotted them in solid black *)

Show[%261,%271,%278,%285]

sol=NDSolve[{
cyc'[t]== ksynth-(a+b*cdk1[t]^n/(ec50^n+cdk1[t]^n))*cyc[t],
cdk1'[t]== (bkgd1+((((cdk1[t])^ncdc25) )/((bigk1^ncdc25)+((cdk1[t])^ncdc25))))kcdc25*(cyc[t]-cdk1[t])-(bkgd2+((bigk2^nwee1)/((bigk2^nwee1)+((cdk1[t])^nwee1)))) kwee1*cdk1[t]+ksynth-(a+b*cdk1[t]^n/(ec50^n+cdk1[t]^n))*cdk1[t],
cyc[0]==45,
cdk1[0]==0},
{cyc,cdk1},
{t,0,tmax}]

(* Plot the solution as a time course *)

Plot[{cdk1[t]/.sol},{t,0,100},Frame->True,PlotStyle->{Darker[Blue]},LabelStyle->{FontFamily->"Arial",FontSize->18},FrameLabel->{Style["Time"],Style["Active Cdk1 (nM)"]},PlotRange->{0,102},AspectRatio->1]

Show[%290,%294,%300,%305,%310,%315,%320,%325,%330,%335,%340,%348,%354,%359,%364,%369,%374,%379,%384,%389,%394]

Export["Time course for Qiong-n=1-June 6-2.pdf",%395]

(* Plot the solution in the phase plane *)

ParametricPlot[{cyc[t],cdk1[t]}/.sol,{t,0,100},Frame->True,PlotStyle->{Black},LabelStyle->{FontFamily->"Helvetica",FontSize->18},AspectRatio->1,FrameLabel->{Style["Total cyclin (nM)"],Style["Active cyclin-Cdk1 (nM)"]},PlotRange->{{0,102},{0,102}}]

Show[%286,%291,%297,%302,%307,%312,%317,%322,%327,%332,%337,%344,%351,%356,%361,%366,%371,%376,%381,%386,%391,%396]